Crystallization of the Lewis- Wahnstrom ortho-terphenyl model 
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Crystallization is observed during long molecular dynamics simulations of bent trimers, a molec- 
ular model proposed by Lewis and Wahnstrom for ortho-terphenyl. In the crystal, the three spheres 
that make up the rigid molecule sit near sites of a body centered cubic lattice, the trimer bond 
angle being almost optimal for this structure. The crystal exhibits orientational disorder with the 
molecules aligned randomly along the three Cartesian axis (an example of cubatic orientational 
order). The rotational and translational mobilities exhibit only modest decreases on crystallization, 
by factors of 10 and 3 respectively. The rotational relaxation does change from Debye-like in the 
liquid to large angle jumps in the crystal. We consider the origin of the superior glass forming ability 
of the the trimer over the monatomic liquid. 



INTRODUCTION 

Molecular crystals, unlike their atomic counterparts, 
must accommodate the enormous variety of particle 
shape, symmetry and flexibility that molecules provide. 
This accommodation typically results in low symmetry 
crystals with a high propensity for polymorphism and 
twinning [I| [2] . There have been relatively few simula- 
tion studies of the crystallization of molecules from their 
melt and, of these, most have focussed on flexible poly- 
mers [3] or high symmetry molecules, either linear [3] 
or octahedral [5]. One notable exception is the freezing 
of water (a molecule with only Civ symmetry) which 
has, understandably, been the subject of a large number 
of studies [5J. With ordering dominated by directional 
hydrogen bonds, water is a special case and, therefore, 
offers little in the way of general insights into the nature 
of ordering of low symmetry organic molecules. 

In this paper we report on molecular dynamics simula- 
tions of the crystallization of a liquid comprised of bent 
trimers (another Civ molecule), inspired by the organic 
glass- former ortho-terphenyl (OTP), in which the inter- 
actions are based on the Lennard- Jones potential. Since 
its introduction in 1955 by Andrews and Ubbelohde [7] 
as a useful liquid for the study of the glass transition, the 
properties of supercooled OTP have been the subject of 
extensive research [8]. In 1993, Lewis and Wahnstrom 
(LW) [9] introduced a model of OTP based on a rigid 
bent trimer of spherical particles which has subsequently 
been used for a number of simulation studies of glassy 
behaviour [TUHI2]- 

Beyond the general similarity of the bent shape, the 
model bears little resemblance to OTP. The trimer an- 
gle, for example, was chosen to be 75° rather than the 60° 
value indicated by the structure of benzene because the 
authors regarded the former value, "a value such that 
crystallization is prohibitively difficult" [10 . The idea 
that it is possible to decide whether a molecule can crys- 



tallize easily or not (however 'easily' might be defined 
here) simply by inspecting the shape of the molecule is, 
we suggest, both widespread and poorly justified. The 
crystallization of the trimer, as reported here, provides 
a salutary illustration of the stubborn ingenuity of pe- 
riodic packing and the questionable value of our current 
intuition about how shape influences crystal stability and 
crystallization kinetics. 

MODEL AND ALGORITHM 

The trimer is a rigid 3-particle complex, the site-site 
interactions all being a spherically symmetric Lennard- 

Jones potential u{r) where u(r) — 4e^(^) 12 — (f) 6 ^ 
with s = k B x 600K ~ 4.988 kJ/mol and a = 0.483 
nm. The two trimer bonds are of length a and the fixed 
bond angle was chosen [TU] to be 75°, presumably to 
avoid the choices 60° or 90°, either of which might have 
been regarded as too easily accommodated within a face 
centered cubic lattice. (Although not really relevant to 
the subject of this paper, we note that real bent trimer 
molecules must have bond angles between ~ 100° and 
120°. Ozone, with a bond angle of 116.78 [13], has been 
proposed by Angell [14] as a potentially good glass for- 
mer based on the low value of its freezing point relative 
to its boiling point.) OTP, in contrast, has a bond angle 
60° and, rather than spherical particles, it is composed 
of benzene rings that resemble oblate ellipsoids. Due 
to steric interactions ('overcrowding') the three benzene 
rings cannot be coplanar and, in the lowest energy con- 
figurations, the three rings lie on different planes. The 
crystal structure of OTP is orthorhombic (i.e different 
unit cell lengths along the three Cartesian axis) with 4 
molecules per unit cell [15], quite different, as we shall 
see, from the crystal structure of the LW trimer. More 
accurate models of OTP have been studied [16] in which 
each of the 18 C-H groups is represented as a spheri- 
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cal particle. These simulations are more time consuming 
than the LW model to run and, to date, no crystallization 
has been reported. 

While the LW trimcr may not provide a particu- 
larly realistic representation of OTP, it is of consid- 
erable interest as a computationally efficient and well 
studied model of a supercooled liquid comprised of a 
low symmetry molecule. We simulated a liquid of LW 
trimers using the GROMACS software package [IT] . The 
tail of the Lennard- Jones interactions were smoothly 
brought to zero at 2.5a using a shift function S: U — 

^ 4e ((4) 12 -te) 6 - 5 (^)) whcrc s ^ = i 

for r < 2.3, S(r) = a(r - 2.3) 3 + P(r - 2.3) 4 + 7 for 
2.3 > r > 2.5 and S(r) = r~ 12 - r~ 6 for r > 2.5 with 
a = 0.2889 . . ., (5 = -0.7789 ... and 7 = -0.005145 . . .. 
Bonds were held rigid using the lincs algorithm [15] , 
Each Lennard- Jones force center were given a mass of 
76.768 au giving a total molecular weight of 230.30 g/mol. 
Trajectories were evaluated for N — 324 molecules in a 
periodic cubic box using a 8 fs time step for the inte- 
grator. Density and temperature were held fixed using 
the Nose-Hoover thermostat with a time constant of 1 ps 
[SI- 
RESULTS 

Spontaneous Crystallization of the Supercooled 
Liquid 

Simulations were performed at several temperatures 
along the zero pressure path, and along four density paths 
(p = {1.066,1.093,1.135,1.149} g/ml) previously inves- 
tigated [T5]. At sufficient supercooling, we observed an 
abrupt decrease in the potential energy from a station- 
ary state (see Fig. [T]) due to crystallization. The average 
lifetime i CI y S t of the supercooled liquid prior to crystal- 
lization at 1.135 g/mL and 375 K was i cr y S t = 2 x 10~ 6 
s, a value roughly 10 2 times that of r Q , the structural 
relaxation time of the supercooled liquid. 

The crystalline order is clearly evident after the drop in 
potential energy as can be seen in the configuration, gen- 
erated at the same conditions as those used in Fig. [T] de- 
picted in Fig. [2} We also observed crystallization at {250 
K and 1.0937 g/mL; 260 K and 1.09272 g/mL; T=275 K 
and 1.09285 g/mL; 375 K and 1.14888 g/mL} (data not 
shown) in simulation runs of ~ 10 _6 s duration. 

Crystal Structure 

Given the low symmetry of the LW trimer, the result- 
ing crystal structure is surprisingly simple. The spherical 
particles that make up the trimer sit near the sites of a 
body centered cubic (BCC) lattice. To understand how 
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FIG. 1: Drop in potential energy due to crystallization of six 
independent runs at p = 1.135 g/mL and T = 375 K. The 
average crystallization time at this density, temperature and 
system size is estimated to t cr y S t ~ 2 x 10~ 6 s. 




FIG. 2: Crystal configuration formed spontaneously from 
the supercooled liquid at p = 1.135 g/mL and T = 375 K. 
The left box display positions of the Lennard- Jones centers, 
that are in a near-BCC lattice. In the right box, only the 
line connecting the two outer atoms of the trimer are shown. 
We note the presence of considerable orientational disorder, 
constrained, however, to orientations that are either parallel 
or perpendicular to each other (i.e. cubatic order). 

this works, consider the idealized arrangements sketched 
in Fig. [3] Note that the smallest possible unit cell of 
the completely filled lattice contains two molecules, an 
example of which is shown in Fig. [4] 

The angle between bonds drawn from two adjacent cor- 
ner sites to the centre is 2 sin - 1 ( 1 /a/3) — 70.5°. This 
would be the lowest energy choice for the bond angle if 
the bond lengths were increased from a to 1.12ct to reflect 
the fact that, in the BCC lattice, they must lie along the 
diagonal of the cubic cell. If we retain the bond length of 
a, then the bond angle that minimizes the lattice energy 
will need to be larger than 70.5°. 

How far off is the LW trimer from the ideal bond angle 
for incorporation into the BCC lattice? To address this 
question we have calculated the approximate potential 
energy per molecule as a function of the trimer angle (p: 
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FIG. 3: On the left, a trimer (indicated in red) with an op- 
timal bond angle of 70.5° occupies the sites of a BCC lattice. 
On the right, the trimer with a bond angle and bond length 
similar to that of the LW trimer is arranged on a near-BCC 
lattice. 




FIG. 4: The unit cell of an orientationally ordered arrange- 
ment of red trimers. Black sites represent the rest of the lat- 
tice sites, which are occupied by periodic images of the trimer 
sites. This forms a herringbone structure on the (Oil) plane. 
Nearest-neighbour interactions are marked. Note that this is 
an idealized structure missing some orientational freedom as 
discussed in the text. 

Consider the orientationally ordered structure shown in 
Fig. [4] which, since it also fully decorates a near-BCC 
lattice, gives a good approximation to the interaction 
lengths in the actual crystal as formed in simulations. 
The cell is constrained to retain the shape shown with 
three cubic subcells, the cubic edge length varies as tp is 
varied, to ensure it corresponds to the distance between 
the end sites of the trimer. 

The interactions of each Lennard- Jones site on the 
trimers with their 8 nearest-neighbours and 6 second- 
nearest neighbours were considered, apart from those de- 
fined by rigid bonds. The 12 distances, r\, represented in 
Fig. [4] by long-dashed blue lines are slightly shorter than 
those represented by the 6 orange dotted lines, The 
2 magenta hashed lines have a fixed length of r$ = a due 
to the imposed requirement to remain cubic. There are 
16 second-nearest interaction distances per two molecules 
(after removing rigid bonds and double counting) , repre- 
sented by solid lines which are the length of the cube 
edge, a, or a little longer for displaced sites. Thus the 
interaction energy per molecule is approximately 

U cry /N ~ 6u(ri) + 3u(r 2 ) + u(r 3 ) + 8u(a) (1) 
where a/a = \Jl — 2cos(y), d/a = cos i? - , n = 
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FIG. 5: . Potential energy of the crystal as estimated by Eq.[T] 
as a function of the trimer angle tp. The shift function S use 
to truncate interactions smoothly in simulations introduce a 
— 72ey ~ 1.8 kj/mol shift of the shown energy. 

y/{%) 2 + d? + {a-d) 2 , and r 2 = ^(f ) 2 + 2(a - d) 2 . 
The estimated crystal energy, shown in Fig. [5j is min- 
imized at tp = 76.4° corresponding to a density of 0.352 
molecules per er 3 or 1.19 g/mL. This suggests that the LW 
trimer angle of 75° is, in fact, a near-ideal geometry to 
crystallize into this dense, low-energy crystal structure. 

While the constitutive atoms are quite comfortably 
accommodated on the BCC lattice, the trimer orienta- 
tion exhibits considerable disorder. The orientation of a 
trimer can be established by the vector parallel to the 
line joining the two outer atoms. In the right panel on 
Fig. [2] we show one configuration of these orientational 
vectors in the crystal (where we have omitted the sense 
of the vector). In contrast to the regular arrangement 
of constituent spheres, the orientational vectors exhibit 
considerable disorder. The disorder is not completely 
random since restrictions are imposed by the BCC or- 
dering of the spheres. 

The result is a form of orientational order known as 
cubatic |20j in which the orientation vectors (as we have 
defined it) are either parallel or perpendicular. This re- 
striction is evident in the peaks in the distribution of an- 
gles between pairs of orientation vectors at 0°, 90° and 
180° shown in Fig. [6] The smaller peaks at ~ 55° and 
125° correspond to orientational defects. The appear- 
ance of this cubatic order during crystallization can be 
monitored through an order parameter defined as 

n 1 ^cos 2 (2%)-^ 



where Oij is the angle between the orientation vectors of 
molecule i and j and is the integral of cos 2 (20 ij) over 
a random distribution of orientations. A global cubatic 
order parameter Q is defined as the average Q,. Note 
that Q is normalized so that zero correspond to random 
orientations, and unity to perfect cubatic ordering. 
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FIG. 6: Distribution of the angles between the orientation 
vectors of neighbors in the crystal (as generated by simula- 
tion) and in the liquid. The green dashed line indicates a 
random distribution. 
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FIG. 7: The drop in potential energy (a) is accompanied by 
cubatic ordering (b). Green crosses indicates cubatic order 
parameter Qi's (Eq. [5| and the black line the global cubatic 
order parameter Q (p — 1.135 g/mL, T = 375 K). Negative 
Qi's in the crystalline phase are associated with orientational 
defects . 



The cubatic order parameters along a crystallizing tra- 
jectory are shown on Fig. [7| We find that the increase in 
the cubatic order follows closely the decrease in the po- 
tential energy of the inherent structures during the course 
of crystallization. In particular, we can see no sign of a 
development of the orientational order preceding the on- 
set of crystallization. 



Relaxation Kinetics in the Liquid and Solid Phases 



In light of the disorder retained by the trimer crystal, 
how mobile are the molecules in the crystal? In Fig. 
[8ja) we plot the autocorrelation function for the angle 
between a molecule's orientation vector at time t = 
and the vector at some later time t for the liquid and the 
crystal phases. We find that rotational relaxation persists 
into the crystal phase and the rotational relaxation time 
increases, on freezing, by only a factor of 10. 

The ratio of the relaxation times of the Legendre poly- 
nomials of rank 1 and 2 is ~ 2 in the liquid and ~ 1 
in the crystal, indicating that the rotational relaxation 
changes on freezing from Debye-like diffusion in the liq- 
uid to relaxation involving large angular jumps in the 
crystal [21] . In Fig. [S^b) we have plotted the distribution 
of angular displacements after a short time (t = 1 ns, 
roughly 5% of the liquid rotational relaxation time) in 
the liquid and crystal phases at density 1.135 g/mL and 
temperature 375 K. While the liquid shows the smooth 
distribution expected for incremental displacements, the 
distribution of angular jumps in the crystal exhibits a 
number of gentle peaks. The predominant peak does not 
occur at an angle of 90° or 180°, as one expect if the 
jumps where from on crystal-aligned position to another, 
but, rather, at an angle of ~ 55°. We conclude that ro- 
tational relaxation in the crystal involves orientational 
defect configurations as transient intermediates. 

The mean square displacement, plotted in Fig.[9j shows 
a similar picture of only a modest slowing down on crys- 
tallization. The centre of mass motion of the trimer in 
the crystal continues to exhibit diffusion-like behaviour 
at long times. At T — 375K the center of mass diffusion 
constants for the liquid and crystal are 1.6 x 10~ 12 m 2 /s 
and 5.3 x 10~ 13 m 2 /s, respectively, corresponding to a 
slow down of only a factor of three. These calculations 
have been performed at a finite fixed volume and there 
will be a pressure drop occurring at crystallization. Even 
taking this into account, the persistent mobility of the 
trimers in the crystal is striking. 

We find that the mobility in the crystal is sensitive to 
the amount of orientational defects present and that, in 
some runs, this defect density decreased through an ag- 
ing process after crystallization. Specifically, in four of 
the crystals formed at p = 1.135 g/mL and T = 375 K 
(see Fig. [I]) the amount of orientational defects fluctu- 
ates around 10%, while in the other two runs, we observe 
aging of the crystal, characterized by an increase in the 
cubatic order parameter. The relaxation times for trans- 
lational and rotational diffusion are observed to increase 
with the increase in orientational order. We emphasize, 
however, that the observed crystallization in all runs in- 
volves a transition to the crystals with ~ 10% orienta- 
tional defects. 
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FIG. 8: (a) Rotational autocorrelation functions defined as 
Ciit) = (P i {cos[6'(0)]}P i {cos[6'(t)]}} where 6 is the angle be- 
tween the trimer base and either x, y or z and Pi is the Leg- 
endre polynomial of either rank 1 = 1 (full lines) or / = 2 
(dashed lines). Ci(t)'s are evaluated for both the liquid (blue) 
and crystal (red) at p = 1.135 g/mL and T = 375 K. Rota- 
tional relaxation time is only affected by a factor of ~ 10 
on crystallizing, (b) Distribution of angular displacements of 
molecular orientations after 1 ns. 
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FIG. 9: Mean squared displacement of LJ centers (full lines) 
and center of masses (dashed lines) in the crystal (red) and 
liquid (blue) phases at p = 1.135 g/mL and T — 375 K. Center 
of mass diffusion constants in liquid and crystal, estimated 
from long time displacements, are 1.6 x KT 12 m 2 /s and 5.3 x 
10"'" 



DISCUSSION 

In this paper we have shown that the Lewis and Wahn- 
strom bent trimer will crystallize. The resulting struc- 
ture is a plastic crystal in which rotational motion and 
translational diffusion persist at rates only modestly re- 
duced from those of the liquid. While the observation of 
crystallization contradicts the model designers' expecta- 
tion that crystallization of the trimer would prove "pro- 
hibitively" difficult, we find that their intuition that the 
trimer is, in some sense, difficult to freeze is still cor- 
rect. To see this, consider the minimum crystallization 
time t™yst f° r the LW trimer and the atomic Lennard- 
Jones liquid. {V^ryst corresponds to the 'nose' in the 
standard time-temperature transformation plot.) For the 
atomic liquid (at a density of p = 1.135 g/mL) we find 
t™yst = 3 x 10~ 10 s while the smallest values of t cryst ob- 
tained so far for the LW trimer is ~ 2 x 10~ 6 s. The factor 
of 10 4 difference in these two values of t™y St corresponds 
to a substantial increase in the glass forming ability of 
the trimer over that of the atomic liquid. 

How we account for this increase represents an impor- 
tant puzzle. It has been suggested that slow crystal- 
lization is a consequence of 'bad' (read 'high enthalpy') 
crystals [55]. The idea is attractive as it reduces the 
problem of adjusting glass forming ability to one of ad- 
justing the lattice energy of the stable crystal. A di- 
rect measure of crystal stability is the change in energy 
on freezing relative to the energy of the metastable liq- 
uid, i.e. AU crys t/U-^Q. At the temperatures associated 
with the respective peak crystallization rates, we find 
AUcryst/Uuq — 0.058 for the trimer and 0.089 for the 
liquid of LJ atoms. 

This preliminary comparison supports the proposition 
that the improved glass forming ability of the trimer 
might be a result of the decreased stability its crystal, 
relative to that of the atomic crystal for which shape is 
not an issue. To establish whether this stability differ- 
ence is sufficient to account for the trimer's glass forming 
ability requires additional analysis that we leave for a 
future paper. 
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Supplementary figure: The Lewis- Wahnstrom ortho- 
terphenyl crystal. 



